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ABSTRACT 



Any change in neutronic properties in a reactor operat- 
ing at steady state will result in a change in the equilib- 
rium neutron flux and hence, the power of the reactor, A 
main cause for a change in neutronic properties is the high 
temperature attained in a reactor, which produces a feedback 
in the reactor operation. The response of the reactor to 
a particular feedback is analyzed by using Liapunov’s 
Second Method to specify stability regimes. Both, the 
point-kinetics model and the distributed parameters system 
are analyzed. 

Data from a typical thermal and fast reactors is used 
to specifically determine the stability domains. 
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LIAPUNOV f S SECOND METHOD 



I . 



A. HISTORY 

In 1892 the Russian mathematician A. M. Liapunov pub- 
lished a long paper dealing with the problem of stability 
of motion ( 1 ) . This paper was translated into French in 
1907 and reprinted in America in 1949. Liapunov 1 s Theory 
received by then only little attention due to the diffi- 
culty in understanding the advanced mathematical theorems, 
to the abstract way it was presented and to its lack of 
practical application and for a long time it was nearly 
forgotten. About 35 years ago, Soviet mathematicians re- 
sumed the investigation of Liapunov’s Theory and its excel- 
lent application in several technical fields, mainly in 
control engineering, was noticed. Significant work in this 
area was published by Malkin ( 2_) , Letov (_3) , Lur 'e (4_) and 
Chetaev (50. The excellent paper by Massera (6^) and the 
translation to English of most of the Russian works stirred 
up considerable curiosity in this country. Bellman (7_) in 
1953 published an excellent work concerning stability, but 
the section on Liapunov’s method is difficult to comprehend. 
According to this author, Hahn (8^), Krasovskii (9_) and 
LaSalle and Lefschetz (10.) have the best treatises on the 
subject, the consolidation of the concepts of Liapunov’s 
Theory and the clear presentation of it, make these three 
books the main references for workers in the field. 
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Recently, Zubov ( 11 ) found the best extension of the 
Liapunov’s Theory, to include the analysis of partial dif- 
ferential equations, which is by now the main topic for 
research. Yet, a good amount of work remains to be done. 
All the previous works, including Liapunov's original 
theory, were devoted to the analysis of stability of 
ordinary differential equations. 

Many authors have applied Zubov's extension to works 
concerning vibrations, reactor physics, hydrodynamics, 
magnetohydrodynamics and control processes. The list of 
references will provide the interested reader with the 
main works in this field. 

B. THEORY 

The name "second method" (or "direct method") is of 
historical origin. Liapunov's Theory is divided in two 
categories. The "first method" which comprises all pro- 
cedures in which the explicit form of the solutions is 
used, specially when represented by infinite series and the 
investigation of stability "in the small" (local stability) 
studies the singular points of a nonlinear differential 
equation by using the appropriate linearized version of the 
differential equation near the singular point. The "second 
method" attempts to make stability statements by using (in 
addition to the differential equations) suitable functions, 
called Liapunov functions, which are defined in the motion 
space. This method which deserves special study due to 
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its inherently advantages over most of the conventional 
methods, investigates the stability n in the large" (global 
stability) . 

Liapunov 1 s Second Method is in essence a more general 
expression of the Hamiltonian (total energy), and it is 
based in the statement that a physical system loses poten- 
tial energy in a neighborhood of a point of stable equilib- 
rium. It is said that it is more general than the total 
energy method, because unlike the energy of a system, the 
Liapunov function, denoted V, is not unique. This is the 
main reason why the "second method" is a tool in the analy- 
sis of stability of dynamical systems, which must be used 
with considerable skill. One of the main features of this 
method is its appeal to geometric intuition; V(x,t) can be 
seen as a measure of the "distance" of the state (x,t) from 
the origin, in the state space and the variation of this 
"distance" (norm of the differential equation) as t varies 
will provide definite bounds of stability regions for the 
prescribed system under consideration. 

Normally stability with respect to one norm does not 
imply stability with respect to another norm. This diffi- 
culty does not appear in finite-dimensional systems since 
all norms defined on a finite-dimensional vector space are ■ 
equivalent . 

1 . Fundamental Concepts and Definitions 

Stability is a property of certain systems of dif- 
ferential equations and is basically concerned with the 
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question of whether or not the dynamic system will return 
to a particular state after it has been disturbed in some 
way . 

Differential equations in their most general way 
can be expressed as 

8y ^ ? t) = F[U(x,t) ,t] 

where U(x,t) is a multidimensional function of space and 
time. It may well happen that F depends upon U(t) alone 
and not upon time explicitly. Then the previous equation 
assumes the simpler form 

U(t) = F [U ( t ) ] 

A system of this nature is known as autonomous. Since it 
is not intended in this work to expose the reader with the 
mathematical aspect of Liapunov’s Theory, in general, the 
presentation of the theory will be made without proofs. 

Let us denote ft(R) the spherical region ||x|| < R and by 
A(R) the sphere boundary ||x|| = R. The matter of concern 
is the stability of the origin. Initiating the motion at 
a point x^, the origin is said to be: 

(a) Stable whenever the path remains in the spherical 
region ft(R); that is, the path never reaches the 
boundary sphere A(R), Figure 1, 

(b) Asymptotically stable whenever it is stable and 
in addition the path tends to the origin as time 
increases indefinitely. 
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Figure 1. Domains of stability. 
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(c) Unstable whenever the path reaches the boundary 
sphere A (R) . 

For autonomous systems: 

A Liapunov function, V(U), is defined as a posi- 
tive definite scalar function with the following properties 

(a) V(U) is continuous with its first partial deriva- 
tives in a certain open region about the origin 

(b) V(o) = 0 only 

(c) Outside the origin and always inside the open 
region, V is non-negative and vanishes only at the 
origin. The origin is an isolated minimum of V; 
and in addition 

V 0 in the open region 

If V = 0, the stability is neutral 
If V < 0, the stability is asymptotic. 

For nonau t onomous systems: 

Let us define V-^(U) as previously, A Liapunov 
function, V(U,t) is defined as a positive definite scalar 
function with the following properties: 

(a) V(U,t) is defined in the open region for all 
t 21 0 . 

(b) V(o,t) = 0 for t -I 0 • 

(c) V^(U) V(U,t) for all x in the open region and 
all t >1 0 ; and in addition 

V <_ 0 in the open region 
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where 



• H V 

V(U,t) = 37- + u • vv 



Summarizing these conditions, the three main 
steps to test if a scalar function is a Liapunov func- 
tion a r e : 

1. V > 0 

2. V (o) » 0 

3. V < 0 
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2 . Methods to Construct a Liapunov Function 

Several authors have attempted to present guide- 
lines for generating Liapunov Functions. Unfortunately, 
most of these methods are extremely restrictive, limita- 
tions being imposed primarily by the number of nonlinear 
terms and the order of the system. Ingwerson (12_) pre- 
sents a Table to generate Liapunov Functions for linear 
autonomous differential equations up to the fourth order. 
Barbashin (1J3) and Lur'e and Rozenvasser (1 ^4 ) tried to 
establish some rules for construction of Liapunov Functions. 
Schultz (1J^) has one of the best methods available now; it 
is called the Variable Gradient Method, being the least 
restrictive for the case of autonomous systems. Its only 
restriction is that all nonlinearities must be single- 
valued . 

In general, literature late in 1961 and 1962 became 
saturated with methods for generating Liapunov Functions, 
most of them with the restrictions previously mentioned. 
Furthermore, such methods, when applicable, were directed 
to autonomous systems only. 

Although development of Liapunov stability theory 
and applications to ordinary differential equations has 
progressed rapidly, its application to partial differential 
equations has remained limited. The importance of partial 
differential equations in the fields of reactor physics, 
hydrodynamics, control processes, etc. has motivated inves- 
tigations of possible ways to extend Liapunov stability 
theory to partial differential equations. 
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In this study Functional Analysis plays an important 
role, because the stability of the equilibrium solution is 
defined in terms of the norm induced by the inner product of 
the Hilbert space on which the solutions of the system are 
defined. Thus, with proper choice of inner product or norm, 
the square of the norm becomes the Liapunov functional which 
establishes asymptotic stability. In fact what is being 
done is the construction of a scalar function of the distance 
between the solution and the equilibrium point of interest. 

If this distance function, evaluated along the solution 
paths, obeys certain inequalities, then statements concern- 
ing the stability of the equilibrium point can be made. 
Movchan (2JL) defines stability in terms of two metrics, 
rather than one, to be more restrictive on the initial 
states. The choice of the initial state space and the 
metric is crucial in the formulation of stability problems 
in the framework of Liapunov’s Second Method. For cases in 
which only the behavior of some of the state variables is 
of importance or some function of the state variables is of 
interest to the analyst, then, for these cases, it is mean- 
ingful to define Liapunov stability with respect to two 
metrics. For instance, one may be only interested in the 
behavior of the maximum deflections in an elastic system, 
regardless of the velocity and acceleration involved into 
attained it . 

According to Kastenberg & Ziskind ( ljO extreme care 
has to be taken when interpreting the obtained results, 
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because severe peaking in some of the state variables of the 
system can cause the norm of the system trajectory to 
move out of the domain* Wang (1_8^ and Parks (19_) have de- 
voted special attention to the study of panel flutter which 
represents a linear distributed parameter system. It is 
concluded that a good comprehension of Hilbert and Banach 
spaces will provide the analyst with an excellent tool for 
generating Liapunov functionals. 

3 . Comparison with Other Methods 

The advantages of Liapunov’s Second Method over the 
conventional methods used for stability analysis can be 
summarized as follows: 

(a) It employs the system equations without resorting 

to approximations. Normally a distributed system is approxi- 
mated by a lumped parameter model having finite or infinite 
number of degrees of freedom. This method is not satisfac- 
tory because quite often it exhibits characteristics which 
do not agree with the physical nature of the problem, 
Liapunov's method is in general a more reliable method. 

(b) There is no theoretical limit on the number of non- 
linearities or on the order of the differential equation to 
which Liapunov's Second Method can be applied. 

(c) The laborious work implied in finding the system 
solution is avoided. Normally the system equations are in- 
tegrated numerically for some given perturbations in the 
initial conditions. This method is time consuming and 
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sometimes does not give an accurate presentation of the 
system behavior. 

(d) The relationship between the system stability and 
the system parameters is directly extracted from the analy- 
sis of Liapunov’s method conditions for stability. 

(e) There are no mathematical restrictions with res- 
pect to uniqueness. 

The method presents also some deficiencies; among 
them the main ones are: 

(a) The construction of the required Liapunov Function 
for the system under consideration is without doubt the 
most difficult part of the task due to the fact that there 
is no guideline to construct it and success depends mainly 
upon the ingenuity and experience of the analyst. 

(b) The interpretation of the obtained results has to 
be done carefully. 

(c) Integral inequalities play an important role in 
deriving conditions for stability. As a consequence, the 
regions of stability may be somewhat looser than those ob- 
tained by other methods. 

In conclusion, in spite of the difficulties cited 
above, the advantages of the Liapunov’s Second Method makes 
it an excellent method for studying stability. 

C. APPLICATION TO NUCLEAR REACTOR CONTROL 

Liapunov’s Second Method has been applied with success 
by Hsu (2_0) who considers the linear and non-linear cases, 
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analyzing them by means of the Spectral Analysis and Liapu- 
nov’s Method, and comparing the results obtained by each 
method. llsu shows that the Spectral Analysis provides a 
little more information than Liapunov's Method, but this 
latter eliminates the calculation of the eigenvalues of the 
system, and for the nonlinear system analysis bounds of 
st-ability are presented. 

Kastenberg (_2_1) presents a clear comparison of Liapunov 1 s 
Second Method with the Semigroup Method, and the Comparison 
Function and Maximum Principle Techniques. In the Liapunov 
and Semigroup methods, one must obtain an a priori bound on 
the system nonlinearity and then proceed. When employing 
the maximum principle or the comparison function technique, 
an a priori bound on the system nonlinearity is not always 
required. In contrast, one must give a bound on the initial 
condition. For cases which are stable in a global sense, 
the results of the various methods coincide. 

The application of Liapunov's Second Method to nuclear 
reactor control seems to be excellent, mainly for the spa- 
t ially-dependent reactor system, due to the fact that the 
method allows the inclusion of any number of nonlinearities. 
This will permit the study of spatial effects of temperature, 
the main feedback effect in a reactor, control rod motion 
and several other processes within a reactor. Also, the 
inclusion in the analysis of the different reactivity coef- 
ficients, as doppler resonance broadening and structural 
expansion is allowed. 
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As pointed before, the interpretation of the results 
must be done very carefully due to the fact that a Liapunov 
Function is not unique. 
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II. MATHEMATICAL MODEL OF THE REACTOR SYSTEM 



Both Thermal and Fast reactors are described essential- 
ly by the same basic dynamic principles, regardless of 
material and geometry considerations. 

For a given reactor size, the reactor reactivity depends 
on the neutron cross-sections and on the relative amounts 
and densities of different materials. All of these being 
affected by the temperature, the reactivity will be strong- 
ly coupled to the power of the reactor and the reactor 
governing equations become nonlinear. 

The Thermal Reactor stability is analyzed by means of 
the lumped-parameter (point-kinetics) model, in which the 
partial differential equations are reduced to ordinary dif- 
ferential equations via spatial discretization. This is 
justified only when infinitesimal small perturbations are 
introduced, as in the case when k is very near to unity 
(very small departures from ’’critical") . Thermal reactors 
are generally smaller in size than fast-breeder reactors. 

The Fast-Breeder Reactor stability is analyzed using the 
governing partial differential equations without resorting 
to any type of approximations. Due to possible stronger 
space effects in fast reactors than in thermal reactors, 
the space and time of the state variables of the system 
must be maintained during the analysis. The reactivity in- 
sertion in the system, Figure 2, can be - either positive or 
negative. Both cases will be considered in this work. 
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A. POINT-KINETICS MODEL 

1 . Diffusion Equat ion 

According to Meghreblian and Holmes (2^2), the neu- 
tron diffusion equation is written in the following way: 



D V 2 f(r,t) - Z a <}>(r,t) = £ - 

- (l-3)v E f pg <J> (r , t ) - Pg^^i C i (r,t) 



( 1 ) 



Expressing the equation for a slab reactor and using only 
one-delayed neutron group: 

D t). _ z d)(x t) » - t - ) - 

dX 



- (l-3)v Z f pg <Kx,t) - pg 1 c(x,t) 



( 2 ) 



It is looked for a solution which is separable in time and 
space : 

<j>(x,t) = <f>(t) F(x) (a) (3) 

C (x , t ) = C(t) F (x) (b) 

Applying relations (3) to Equation (2) yields 



F(x) 3x 2 



9 2 F(x) _1_ 

a <|>(t) 



- <Kt) - 



- (i-3)v E pg 4> ( t ) - pg X c(t) 



Thus 



1 5 2 F (x) _ 1 



F (x) 



3x^ 



-5 K + W) [v i(c) - 



- (l-6)v pg 4>(t) - pg X c(t) J | = - 



(4) 



(5) 
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From Equation (5) it is obtained that: 

+ b 2 f(x) - o 

9x 



(a) 



( 6 ) 



(db 2 +E )4>(t) = - ^4>(t) + (i-3)vZ f pg^(t)+ P g X c ( t ) (b) 

a V I 



Using the following relationships: 



L 2 - 



D 



k = 



V E f pg 

2 2 *™ 
Z (1+L B ) 
a 



% = 



v(Z +DB ) 
a 



and k = Ak + 1> Equation (6b) becomes 



*4>(t)- [ (l-P)Ak-ei(Kt) + — 2^-2 x C ( t ) (7) 

E +DB 
a 

Let P(t) be the time behavior of the reactor power 
generation 

P(t) = e Z f (Kt) (8) 

• • 

then P(t) = e Z^ (j)(t), and the equation describing the 
reactor power is 

£ £ f pg 

£P(t)= [ (1-g) Ak-3 ] p ( t) + 2 * CO) (9) 

Z + DB 
a 

It will be seen in the next paragraph that for certain 
kinetics problems a simplification can be obtained by 
assuming the infinite delay time in delayed neutron emis- 
sion: 

XC(o) = P(o) 
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This expression is used as a first estimate for the precur- 
sor concentration. 

Then 



£P(t) = [(l-B)Ak-B] P(t) + [ 1+Ak ] B P(o) 



( 10 ) 



2 . Concentration of Precursors Equation 

yjr C i (r,t) = - \ ± c j ,(r,t) + vB ± 2 f <j)(r,t) (11) 

For a slab reactor and using one-delayed neutron group, 
Equation (11) can be expressed as 

•|y c(x,t) = - X C (x , t ) + v B £ f <f> (x,t) (12) 

Using Equations (3) to separate variables, Equation (12) 
becomes 

C ( t ) = - X C(t) + V 3 2 f <KO (13) 

Then, using Equation (8) yields 



c(t) = - X c ( t ) + ^ B P(t) 



(14) 



From Equation (14), the steady state concentration of the 
precursors can be obtained as: 

X C (o ) = P (o) (15) 



This result was previously used. 
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3 . Energy Equation 



Expressing H(r,t) as the energy content per unit 

volume, then the time rate of change of H must be given 

by the net energy gain per unit time and volume from the 

fission reactions and the reactor cooling: 

H ( r , t ) = pC [ T (r , t ) - T (r) ] = pC 0(r,t) (16) 

P op 

and 

— - H ( r , t ) = £ E f <j> ( r , t ) - (P ( r,t) (17) 

For a slab reactor, the energy balance becomes 



PCp |^- 6 (x, t ) = e £ f e (x, t) -/^(x,t) (18) 

Let 0(x,t) = 0(t) F (x) (a) (19) 

4>(x,t) = <KO F(x) (b) 

(? (x, t) = (p (t) F (x) (c) 

Then 

pc p 0 (t) = e £ f 4>(t) - (t) (20) 



Let C f 



pC and using Equation (8), Equation (20) becomes 



C 0(t) = P(t) - (?(t) 



( 21 ) 



Three cases are normally encountered for the power removal: 

(a) Adiabatic, CP (t) = 0 

(b) Constant power removal, (p (t) = A 

(c) Newton’s law of cooling, (p ( t) = h0(t) 
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Feedback Models 



4. 

a . Temperature Feedback Model #1 

The multiplication factor is temperature and 
time dependent 

k = k(T) = k [ T ( t ) ] (22) 

Expanding it in powers of [T(t) - yields 



t(I > -e 



oo n 

dk(T) 



dt 1T(t) - V 



n 



n = 0 



thus 



E 

n = 0 



k(T) = k Q [1 + >: An (T T q ) ^ 



n 



(23) 



which becomes, upon dropping all terms beyond n = 1, 
k = k Q [l + Y (T~T o ) ] = (1 + Ak Q )(l + y0) 
or 




(24) 



(25) 



where Ak^ is the change in k applied to the reactor at 
t = 0 and y represents the temperature coefficient of 
reactivity, which could be either (+) or (-) . 
b. Feedback Model #2 

One of the major considerations to be made when 
accounting for the temperature dependence of the multipli- 
cation is that described by the nuclear Doppler effect. 

Thompson and Beckerley (2J^) expresses that 



dependence as 



dk _ . 300 , n 

dT ~ Y D^ T ; 



(26) 



29 



where T is in °K. 



From (26), Ak can be expressed explicitly as 



Ak = 



Ak + 
o 



Id_ 

1-n 



ETC 



300 , n 
T 



- V 



300 v n 
T 

1 o 



(27) 



being n = 1 a special case, for which a logarithmic expres- 
sion is found. 

Typical values of n are 1/2, 1 and 3/2. 



B. DISTRIBUTED-PARAMETER REACTOR SYSTEM 
1 . Diffusion Equation 

The diffusion equation is used in this case without 

resorting to any approximation, then 

DV 2 <J>(r,t) - Z a <J>(r,t) “ ^ f^r (r,t) - 

-(l-£)vE f p g <J>(r,t) - PS d C ± (r , t ) (28) 

i 

For a slab reactor and using one-group delayed neutron, the 

diffusion equation takes the following form: 

E— |(x,t) - <f>(x,t) = ~ |^-(x,t) - 

3x 

-(l-3)v E f pg <{>(x,t) - pg X C (x , t ) (29) 

The cross sections should be written in a proper way as 

time-dependent but these variations are assumed negligible 

in comparison to the rapid transients in <f> , C and T. The 

equation for a non-trivial, steady state solution cf) ^ ( x ) » 

C (x) is : 
o 

3 2 <j> (x) 

D x - E d> (x) = -(l-B)vE f pg <J) (x)-pg X C (x) (30) 

3x a ° to o 
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Now subtracting Equation (30) from Equation (29) yields 



D— f (x, t) 
3x 


- Z a \p(x, t) = ~ If (x, t) - 




- (1-8) v Z f 


pg ^ (x , t ) - pg X -f (x,t) 


(31) 


where 


ip(x,t) = <{-(x,t) - <f> Q (x) 


(a) (32) 




(x, t) = C (x, t) - C Q (x) 


(b) 



This change of variables moves the equilibrium state of the 

system from (<f> , C ) to (0, 0) , Figure 3. 

o o 

V E f pg 

Define k = ~ (33) 

I a., 
a 1 

2 2 

and using a = 1 + L B (34) 



the diffusion equation becomes 



_ I ^(x,t) + (l-e)k Z a Tp(x,t) + 



3x‘ 



+ ps X iS(x,t) - 



(35) 



v 3 1 

Dimensionless space and time expressions are introduced at 
this p o in t : 

(x > > L) 



x 

y = _ 



T = vE t 
a 



(a) (36) 

(b) 



Then Equation (35) becomes 




(37) 



The boundary conditions are: (a) ^(W.t) = 0 

(b)tK-W,T) = 0 



x 



(38) 



where W = expresses the dimensionless half-width of the 
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Figure 3. Translation of equilibrium states. 
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slab reactor including the extrapolation distance. The 
initial conditions are: (a) ^(y,o) = 0 (39) 



(b) (y , o) = 0 

2 . Concentration of precursors equation 

3C.(r,t) 

— ^ = -X i C..(r,t) + v3 i E f <J) (r , t) (40) 

For a slab reactor and using one-group delayed neutron, 
Equation (40) is written in the following way: 



9 C (x , t ) 
3 t 



AC.(x,t) + v3 Z^ (|)(x,t) 



(41) 



The same assumption used in the diffusion equation, with 
respect to the cross sections is used here, A non trivial, 
steady state solution <j> o (x), C Q (x) is assumed: 

0 = - XC^(x) + Vg Z^ 4> (x ) (42) 

Subtracting Equation (42) from Equation (41) yields 



9 $ (x, t) 

at 



(x,t) + Vg E f \p(x,t) 



Using Equations (36), Equation (43) becomes 



(43) 



vZ 

a 



9^ (y,T) 

dr 



( y,T) + vg E f \p(y,r) 



(44) 



The initial condition for the equation is 

T^(y,o) = 0 (45) 

Introducing the multiplication factor in Equation (44) yields 



vE 



9^ (y ,t) 

9t 



gaE k 

^ (y»i) + - p - g" ■ 'l'(y.T) 



( 46 ) 
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* 



Figure 4. Dimensional 

slab reactor . 



Figure 5. Dimensionless 
slab reactor . 
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3 . Energy Equation 



The energy balance condition is normally expressed 



as : 



Energy produced - Energy removed 



Energy stored in 
the system 



The energy produced is due to fission: e E^ <f>(r,t) 

The energy removed is due to the coolant: (p (r,t) 

Three cases are normally encountered for the energy removed 
from a system: 



It is noticed that these cases are for stationary-fuel 
reactors. The energy equation can be written as: 



expressed as energy per unit volume per unit temperature. 
The subscripts f and C stand for fuel and coolant respec- 
tively. Using the same assumption as used previously with 
respect to the cross section, and expressing Equation (47) 
for a slab reactor, yields 



(a) Adiabatic case, CP (r,t) = 0 



(b) Constant power removal, for which the steady 
state requirements are normally used, CP (r,o)= eE^(j)(r,o) 

or (p (r) = eE f <j> o (r) 



(c) Newton's law of cooling, CP (r,t) = h 0(r,t) 




3T f (r,t) 



(47) 




C f — 9t = e £ f 4>(x,t) - (p (x,t) 



3T f (x, t) 



( 48 ) 
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Adiabatic Case 



3T (x,t) 

C f — = e E f 4>(x, t) (49) 

This case will not be considered in this work because it 
represents accident conditions* 

b. Constant Power Removal 

3T f (x,t) „ 

Cf = e 4>(x,t) - 0 J Q (50) 



Assuming a steady-state, nontrivial solution ^^(x) it 
becomes 

0 = e E f 4>(x,t) - (P (51) 

Subtracting Equation (51) from Equation (50) yields 

Cf = e Zf ip(x,t) (52) 

where 9(x,t) = T^(x,t) - T^ q (x) (53) 

Introducing Equations (36), Equation (52) becomes 

C f vZ a fr ^ y,T ^ = c TpCy,?) (54) 



with initial condition: 



o(y»o) = o 



(55) 



c. Newton's Law of Cooling 
3T (x,t) 

Cf = c If <j>(x,t) - h(T f - T c ) (56) 

Assuming a steady-state, nontrivial solution, <f> 0 ( x ) and 
T^ q (x), Equation (56) becomes 



39 (x , t ) 
f 3 1 



e E f ^ (x , t) 



h9 (x , t ) 



( 57 ) 
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Note that and h remain constants. This is explained by 

means of the one effective coolant temperature assumption. 
Introducing the dimensionless expressions for space and time. 
Equation (62) becomes: 



c f v z a " e z f - hG (y» T ) 

with initial condition: 0(y,o) = 0 

4 . Feedback Models 

a. Temperature Feedback Model //I 

Through the temperature, the multiplication 
factor is now dependent on both space and time. 



(58) 

(59) 



k = k(T) = k [ T (x , t ) ] 



(60) 



Expanding it in powers of [T(x,t) - T^(x)] yields: 



0 ° 

k(T) = - 3 - - H IT ) [T(x>t) _ t (x ) ] n 

n 9T 
n=0 



Thus 



k ( T ) = k Q [l +£ A n (T - T Q ) n ] 
n = l 



(61) 



which upon dropping all terms beyond n=l, reduces to: 



k = k Q (l + y(T - T q )] = k Q (l + y6) 



where 



Y « 9k/ 9 T = A: 



(62) 



and 0 (x , t ) = T ( x , t ) - T (x) 

o 



Then, Equation (62) becomes 



k = 1 + Ak + Y® 
o 



( 63 ) 
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having neglected the term yAk^0 . The multiplication factor 
can now be written as: 



k[0(y,T)] = 1 + Ak o + y0 (y , T ) (64) 

y is the temperature coefficient of reactivity, which could 
be either (+) or (-) . 

b. Feedback Model #2 

The temperature dependence of the multiplication 
factor can be expressed in a general form, adding the con- 
tributions due to the Doppler effect, due to structural ex- 
pansions and due to other effects, as: 

dk , l,n . ,T . m , / c <z\ 

« ■ Vt“> + + i (65) 



where D and E stands for doppler and expansion respectively, 

k can be expressed explicitly, using the initial condition, 

k = k at T = T , thus 
o o 



k - 1+Ak + 



Yta a i a i 

, r T (— ) n - t (— ) n l + 

o ( 1-n) 1 4 ; o4 ; J 

o 



4 >“ - ? + - T ) 

1-rm a 2 o a 2 o o 



( 66 ) 



when n = 1, the third term on the right hand side becomes: 

(66a) 



Y b a 1 ■ in f 



It is noted that a^ = - 300°K and T is usually expressed 

in °K. 
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III. STABILITY ANALYSIS 



A. POINT-KINETICS MODEL 

The reactor under consideration for this analysis has 
the following characteristics: 

(a) Thermal reactor. 

(b) Homogeneous, bare, slab reactor. 

(c) One-group delayed neutron. 

(d) Constant power removal. 

(e) Stationary - full reactor. 

and normal operating conditions (no accidents) are assumed. 
The governing equations for this system are: 

£P(t) = [ (1-3) Ak -3] P(t) + (1+Ak) 3 P(o) (10) 

C0(t) = P(t) - t) (21) 



using the following feedback models: 



(a) 


Ak = Ak 

o 


+ Y 0 


(b) 


Ak = Ak 

o 


+ 1 



IT<^V 



I <^V) 



(25) 

(27) 



o 

where Ak^ is a positive step insertion of reactivity by 
external means, such as control rod motion. The initial 
conditions for the system are: 



(a) P (o) = (P(o ) 

(b) 0 (o) = 0 

(c) 0(o) = 0, T ( o) = T q 

For the constant power removal, (p (t) becomes P(o) and 
Equation (21) can be written as 

c0(t) = P(t) - P(o) (21a) 
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Assuming P(o) = 1 and introducing a dimensionless power, 
Equations (10) and (21a) can be expressed as 



JlP(t) = [ ( 1 - 3 ) Ak -B]P(t) + ( 1 +Ak) 3 (10a) 

c0(t) = P(t) - 1 (21b) 



1 . No Delayed Neutrons, No Reactivity Input 

The point kinetics model when 8 = 0 and Ak = 0 is 

o 

reduced to 



JiP = yQP (67) 

c0 - P - 1 ( 68 ) 

After some algebraic manipulations, Equation (67) can be 
expressed as 



d 

dt 



In 



P 




(69) 



"Cross multiplying" Equations ( 68 ) and (69) yields 

• d Y c . 

P - 3 — In P - 7 — 00 = 0 
dt i 



(70) 



which can be written as 



f^r - In P - i e 2 ) - 0 



(71) 



The expression inside the parentheses can be called V, 



V = P - InP - 



1 

2 



Y 



^ 0 2 



then 

(72) 



and therefore -r — = 0 . 

d t 

For V to be a Liapunov function, it has to fulfill the fol- 
lowing conditions: 



(a) V (P , 6 ) > 0 for P>0 and 0>O 

(b) V ( 1 , o ) = 0 (73) 

(c) V(P,0) < 0 for P>0 and 0>O 
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If V of Equation (72) is modified as 



1 Y c 2 

V = (P - In P-1) - j e 



( 74 ) 



it is seen that Equation (71) is not altered, and for a 
negative reactivity input Equation (74) becomes 



The expression (P - In - P-1) > 0 for all values of P > 0. 

Then V satisfies condition (73a) for all values of P and 
0 > 0, when y is negative; for a positive value of y, con- 
dition (73a) is satisfied whenever 



Condition (73b) is clearly satisfied for the steady state 
solution P = 1 and 0=0. Condition (73c) gives V = 0 for 
both cases of positive and negative y, which means that 
asymptotic stability has been ruled out. A similar result 
was found by Ergen and Weinberg (2_4) using the Hamiltonian 
approach. Due to its simplicity, this case has been added 
to this work to acquaint the reader with a general view of 
the necessary steps in the Liapunov’s Method formulation. 

This case does not represent any meaningful practical 

# ' > 

s x t uat ion . 

Typical data for a thermal reactor is obtained from 

Solomon and Kastenberg (2_6^) , as shown in Table 1. From 

this data, £ is found to be 1.235 x 10 ^sec. An average 

3 

heat capacity for liquid light water is used, 1 cal/cm - °C. 




(75) 



1 Y c 2 

(P - In - P-1) > X — - 0 
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TABLE 1 






V l r 


0.49780 


-1 

cm 




f 

Z 


0.261 


-1 

cm 




a 

V 


2.21 x 10 5 


-1 

cm-se c 




6 


0 . 0064 


- 




D 


0 .409718 


cm 




a 


1.41 


- 




A 


0.080 


-1 

sec 
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The case of* negative y has been found to be stable 
for all values of P and 0, but the stability is only 
neutral (V = 0) . 

The case of positive y is represented in Figure 6, 
in which the system is seen to be unstable for most of the 
operating ranges of P and 0, 

2 . One-Group Delayed Neutron, No Reactivity Input 
In this case 3 is included in the analysis and 



Ak^ = 0, the governing equations become 

£P = [ (1-3) Ak-g ]P + (1+Ak) 3 (76) 

c0 - P - 1 (77) 

Ak = y0 (78) 

Substitution of Equation (78) in Equation (77) yields 

£P = [(1-3)Y0-&]P + (l+y0)3 (79) 



"Cross multiplying" Equation (79) and Equation (77) gives, 



after some algebraic manipulations 

P- In P = I 1 (1-8)00- ^0+| (1+Y0) — ~ ~ (80) 

Thus 

— -[P-InP- | |(1-B)y0 2 + ^ 0] = | (1+Y0) ■ (P j 1) (81) 

Then 

V = (P - In P-1) - \ ft ( 1 -^)9 2 + y- 0 (82) 

V = f (1+Y0) - ( - P - p 1 -- (83) 



It can be seen that the inclusion of the delayed neutrons 
complicates the Liapunov function, obtained by the same 
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p 




Figure 6. Stability domains for the case of no delayed 
neutrons, no Ak input, positive temperature coefficient 
of reactivity* 
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steps as for the previous case, but this case approaches a 
more real situation. 

Two cases are analyzed: 

a. Negative reactivity: y = - |y| 

V = (P - In P-1) + J (1-3) 0 2 + e (84) 

V = | ( 1- i Y | Q ) (85) 

It is seen from Equation (84) that V > 0 for all positive 

values of P and 0. Normally (l-|y|0) > 0, then for values 

* 

of 0 c P 1, V is less or equal than zero. Clearly V(1.0) 

= 0, then V only vanishes at the equilibrium state. It is 

concluded that asymptotic stability is obtained whenever P 

be less than 1.0, and for all 0 > 0. Since Ak = 0, P can- 

o 

not be larger than 1.0. Figure 7 shows the domain of 
stability. 

b. - Positive reactivity: y = |y| 

Equations (82) and (83) represent this case. The condition 
of asymptotic stability is obtained whenever 

3 i 7 

(P - In P-1) + 0 > j f (1-8)Y0 

and P < 1. 

This case presents a very narrow region of stability, as 
seen in Figure 8. 

3 . One-Group Delayed Neutron, Step Reactivity Input 

The two previous cases have only theoretical inter- 
est because they do not represent a practical situation. 

The case of Ak^ input represents a real practical situation, 
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p 




Figure 7. Stability domains for the case of delayed 
neutrons, no Ak input, negative temperature coeffi- 
cient of reactivity. 



P 




Figure 8, Stability domains for the case of delayed 
neutrons, no Ak input, positive temperature coeffi- 
cient of reactivity. 



4 6 



i.e. start-up and shut-down, planned or unplanned motion of 
control rods. The governing equations for this system are: 



ZV « [(l-3)(Ak + y0)-3] P 4- (14-Ak Q + y0)3 (86) 

c0 = P - 1 (87) 

Substituting Equation (87) in Equation (86) yields, after 
simplifying terms: 

P - [Ak o (l-3)-3] f 0 = i (Ak o + 3y9) + 'sT (1-3)p (88) 

Then 

[Ak Q (i-3)-3] J e] = |(Ak o +3ye)+ J(i-3)0p (89) 

Thus 

v = (P-1) + [3- Ak o (l-3) ] f 0 (90) 

v = j(Ak Q + 3y0) + \ (1-3) 6P (91) 

For negative y, Equation (91) remains unchanged and Equation 
(90) becomes 

V = f(Ak o -3 | Y | 9) - -bd-(l-3)0P (92) 



In order to obtain a positive definite V, it is required 
that 

P + [3 -Ak Q (1-3) ] f 9 > 1 (93) 

This condition is satisfied whenever 

3 



Ak < 



(94) 



“o 1-3 

for all 0 < 0 , in agreement with the linear theory* To get 
a negative definite V, it is required that 

Ak o < [3 + (1-3) p ] I Y I 0 < 0 (95) 



This is the additional stability condition required by the 



47 



nonlinear theory. Equations (94) and (95) are specifying the 
bounds for asymptotic stability for a negative y. It is 
noted that immediately after the step change in the multipli- 
cation factor, the neutron flux increases rapidly; this is 
referred to as the prompt jump, Lamarsh (25), and it is ex- 
pressed as 

3 (1-Ak ) 

p(t) - — g^Air~ p(o) (96) 

p o 

For positive y, it is seen from Equation (91) that V will 

always be positive definite for all values of P and 0, 

representing an unstable situation. Data from Table 1 is 

used to specify domains of stability for various Ak^ inputs. 

The first bound for all these cases is provided by Equation 

(94) giving a Ak Q < 1.00625$, in order to obtain a positive 

definite V, for all P > 1 and 0 > 0. It is clearly seen 

that' P has to be greater than one, the steady state power, 

after a Ak insertion. The second bound is obtained, set- 
o 

ting V = 0 in Equation (92), Figures 9 and 10. After the 
insertion, the power will increase to an amount specified 
by Equation (96). From Figures 9 and 10, it can be seen 
that right after the Ak^ insertion, the system is unstable, 
and P and 0 increase, until the boundary V=0 is crossed, 
then the system becomes asymptotically stable and the tra- 
jectory returns to the equilibrium position, but when V=0 
is crossed again, the system turns to be unstable, and in 
that way, the trajectory keeps oscillating along the line 
V=0. A similar situation is obtained by Meghreblian and 



48 




Figure 9. Stability domains for the case of delayed 

neutrons, step Ak = 0.156$, negative temperature 

o 

coefficient of reactivity. 
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Figure 10. Stability domains for the case of delayed 
neutrons, step Ak = 0.3125$, negative temperature 
coefficient of reactivity. 
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Holmes (2_2) when solving the system of differential equa- 
tions. Their results showed that this is the case of damped 
oscillation, the presence of the delayed neutrons providing 
the damping force, thus, following an insertion of reacti- 
vity, the system will reach a new steady state power. 

4 . One-Group Delayed Neutron, Ramp Reactivity Input 



where a represents the rate of insertion of reactivity. 
Combining Equations (86) and (87) yields, after some alge- 
braic manipulations 



The governing equations for the system are Equations 



(86) and (87), with the difference that now Ak^ is a func- 
tion of t expressed as 



Ak (t) = at 
o 



(97) 



JLP -((l-B)yc 00+ 3c0 -(l-B)atP - Bat = y6 



(98) 



Then 




(99) 



From Equation (99), V and V are obtained to be: 




( 100 ) 



V - (1-6) P + p- + 



( 101 ) 



For a negative y, Equations (100) and (101) become 



V = P 




( 102 ) 



v = (1-6) y- p - p- + 



(103) 
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The following condition is obtained for asymptotic 
stability : 

Ak o (t) = at < -(Jlfyp+3 (104) 

It can be seen that V is positive definite for all P > 1 
and 0 > 0. Equation (104) is necessary to obtain a negative 
definite V to ensure asymptotic stability. For a positive 

Y, the system is inherently unstable. Equation (101) gives 

♦ 

a positive definite V for all values of P and 0 for t >. 0 . 

A typical ramp insertion of reactivity is shown in Figure A. 
During the reactivity insertion, the system is unstable, 
then V = 0 is obtained after the ramp insertion has ended 
(10 sec) and the power and temperature have risen to an 
appropriate level. Setting V = 0, in Equation (103), the 
following expression is obtained 



Me - BAk o 

(1-3) Ak 

O 



(103a) 



At the end of the Ak insertion (t « 10 sec), this expres- 

o 

sion implies that P depends linearly on 0, Figure 11. 

5 . One-Group Delayed Neutron, General Reactivity Input 
The equations for this system are: 

IV = [ (1-B) Ak-$ ]P + (1+Ak) 3 (106) 

c0 = P - 1 (107) 

Y n 



Ak = Ak [ T (— £— ) n - T (^) n ] 



(108) 



Let 



R(T) = CT(^p) n " T o (^p) n ] 

O 



(109) 
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k 




— > 



t 



Ak = O.OOOlt 
o 



Figure A 
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Figure 11. Stability domain at the end of a ramp 
insertion of reactivity. 
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Then Equation (106) becomes 

ZV = [ (1-6) (Ak o +y D -R(T)-B]P + [ 1 +Y d • R(T) ] 3 (110) 

Combining Equations (110) and (107) yields, after some 

algebraic manipulations 

• • • 

Z? - (1-3) Ak c 0+3 c 0 = [ (1-3)P + B]Y TA R+Ak (111) 

o Do 



Then 



d 

d t 



e-d-6)Ak 

[P - 1 + £ c0] = 

Y n 'R Ak 

[ (1-3) P +3] -j— + -f- 



( 112 ) 



Therefore 

V = P-1 + [3 - (1-3) Ak ] 

o /6 



(113) 



and 

y *R Ak 

V = [ (1-3) P +3] -j + 



(114) 



Equation (113) is independent of y; therefore, V is positive 
definite whenever conditions (93) and (94) are satisfied. 

The similarity with the step insertion of reactivity can be 
seen. Equation (114) becomes, for negative y 



V 



Ak 

o 

Z 



Iy d I *R(T) 

[3 + (1-3)P] — — 0 



(115) 



Equation (115) is negative definite whenever 

Ak o < [3 + (1-3)P] I Y D I - R (T) (116) 

for all n ^ 1. The case n = 1 is the logarithmic case and 
Equation (116) becomes 



Ak Q < 300 [ 3+ (1- 3) P ] 1 Y d | In | 



(117) 



The case of positive y leads to an unstable situation due 
to the fact that Equation (114) is positive definite for 
all n, P and T. Data for this case are found in Thompson 
and Beckerley (2J3) and are typical of a fast reactor. The 
following specific case is analyzed: 

(a) Oxide reactor with volume ratio U0 n to P 0 n = 7 

2 u 2 

(b) Sodium density = 50% 

(c) Y d = 1.28 x 10“ 5 /°K 

(d) n = 0.96 

(e) 3 = 0.0033 

(f) Ak = 0.156$ 

o 

The domain of stability is shown in Figure 12, where the 
V = 0 is obtained from Equation (115). 

B. DI S TRIBUTE D- PARAMETER REACTOR SYSTEM 

The reactor under consideration for this analysis has 
the following characteristics: 

(a) fast reactor 

(b) homogeneous, bare, slab reactor 

(c) one-group delayed neutron 

(d) Newton f s law of cooling 

(e) stationary-fuel reactor 

and normal operating conditions (no accidents) are assumed. 
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340 



Figure 12. Stability domains for the case of delayed 
neutrons and general reactivity insertion. 
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The governing equations are: 



8 _ ^ _ (^T _ ) _ +[(l-3)ka-l]fl(y,T) + A^(y,x) 

3y 2 VZ f 

: 3^(y>x) 



vX 

a 



ijj (v.t) 

8t 



ga £ k 

-X^(y,T) + — tKy.T) 

r c» 



( 37 ) 

(46) 



c f vE a ~ °5 t' , t) = e E f - h0 (y. T ) 



(57) 



using the following feedback models: 

(a) k ( y , t ) = 1 + Ak Q + y0(y» T ) 

y a a 

(b) k = l+Ak + y— — [ T (y~) n - T (—•) n 3 + 

o 1-n 1 o I 

o 



(63) 



~- [T( ~-) ra - T o (^) m ) + Y (T - T ) 
1+m a 2 o a 2 o o 



( 66 ) 



where Ak^ is an external positive reactivity insertion. The 
boundary conditions for the system are: 

(a) ip(± W, t) = 0 
The initial conditions are: 

(a) iKy.o) = 0 

(b) £ (y , o) = 0 

(c) 0(y,o) = 0 

The case n = 1 will give the logarithmic term in the doppler 
expression . 

1 . One-Group Delayed Neutron, Step Reactivity Input 

To formulate the concept of stability a distance be- 
tween the equilibrium position (cj) , C , T ) and any perturbed 
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position (<J), C, T) is defined. The variation of this norm, 
which is induced by the inner product of the Hilbert space 
in which the solutions of the system are defined, with time 
will provide information concerning the perturbed position. 
Let W 



[ d C4) , 4> o ; C , C o ; T, T o > ] 2 = j (A X V+ 2 ) dy 



where 



-W 



p 2 v 2 
2 e 



and 



e 2 v 2 £ 2 



23 2 v 2 



a 2 r 2 2 V 2 
A- = C, v £ 

3 fa 



t = <P - <P q 

= c - C 
r o 

0 = T - T 



The Liapunov function is then defined as 

VO, £ , 6) = | | d | | 2 = < d , d > 



(118) 



(119) 

( 120 ) 

( 121 ) 



( 122 ) 



which by definition is positive definite. It is noticed 
that the introduction of the constants A^ , and A^ makes 

V represent the total rate of energy increase per unit 
volume of the system, which is defined as the distance be- 
tween the equilibrium and the perturbed states. 

V vanishes only at- the equilibrium position 
T q ) ; therefore, the requirements of a negative definite V 
will provide the domain of asymptotic stability. 
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Thus 



W 



dV 

dT 



C£ 2 X f 2 # + ¥ + 2C f Z v Z Z Z e m dy 

f ° y 3 2 v 2 r 3y f « 3v 



2 2 v 2. 86. 
- c v E 0 - 5 —. 
f a 8 y 



-W (123) 

Substituting Equations (37), (46) and (57) in Equation (123) 

yields 



W 



dV 2 2 

d T " E 



-W 



3y 



+ [ (l- 6 )ka-l)tp 2 + (. jp + 



vE 2 ka 



v£ X n 2C^E ka 2C^v£ h 0 

a — r <e + . e» f— ; *- e 2 dy] 



Bv 2 z f 2 



£ 2 v 2 E f 2 



veE f 2 



e 2 E ^ 2 



(124) 



Integration by parts of the first term and using the 

boundary condition \j) (+ W, t) = 0 yields 
W W 

i> dy = - I (|^) 2 dy (125) 

-w -w 




Knops and Wilkes ( 2_7_) , in their work, provide the inequali- 
ties useful for this kind of analysis. The so-called 
"eigenvalue inequality 11 plays an important role in this 
specific problem. For a system defined by the eigenvalue 
problem 



V 2 y + Xy = 0 

in domain with y(a,t) = y(-a,t) = 0 , 
following inequality can be stated: 

a a 

x / v 2 By < | <f ^> 2 

-a -a 



on the boundary, the 



dy 



( 126 ) 
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an inequality that can be proved using the calculus of 
variations. (See Appendix A). In applying this inequality 
to this work, it is assumed that the perturbed reactor has 
the fundamental eigenvalue not appreciably different from 
that of the unperturbed reactor, then 



7T 



2 




X 



0 



(127) 



4W" 



Thus, Equation (126) becomes 



W 



W 




'P 2 dy 



(128) 



Substitution of this result in Equation (124) yields 




W 



2r vE 2 ka 
f a 



2C^v£ h 



f_i!_ 02 



dy 



(129) 



veE f 2 
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. 






After introduction of the feedback term, the follov:ing 
expression is finally obtained: 



,W 



dV 2 v 2 

3 ? i - e z t a 



-w 



[1— (1 — 3) (l+Ak Q ) ] \p' 



y vZ 

+ ( 1 - 3 ) ye^ - (i+Ak )(-f- + — 2 

° vE f 3v 2 £ f 2 v 



‘V X! ^ 

+ (—4— + ) y9 

vE f 3v 2 E f 2 ^ 



vE X 
a 



$ v 2 E 2 a 

f 



2C,vI h 2C f vE 

+ — 6 2 — (1+Ak ) 8 \p + 

. 2 v 2 „ v 2 , . O 



e E f a 



eE f v 



2C ,vE 

+ — f — a - Y e 2 ^ 

eE 2 v 



dy 



Let 



a = 1 - (1-3) (1+Ak ) 
1 o 



vZ X 

a 



32 3 2 v 2 E f 2 a 



2C c vE h 
f a 

^ o 

3 e 2 E f 2 a 



^ vZ 
A , a 

a , = — rr- + 



4 vE f 3v 2 E f 2 



a 5 = 



2C,vE 
f a 

eE 2 v 



(130) 

(131) 

(132) 

(133) 

(134) 

(135) 
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The Equation (130) can be expressed as 

W 



dV 2^2 

— < - e a 

dx — f 



j [ 3-1 <P 2 + Z 2 ft 



2 + a 3 02 + 



-W 



+ (1-3) y0^ 2 - (l+Ak^) a^ + a-^ y0i£ijj - 



- (l+Ak Q ) a 5 Q\p + a 5 y0 2 ^] dy 



( 136 ) 



To obtain asymptotic stability, it is necessary to prove 
that the expression inside the curly brackets is positive 
definite. A strategy successfully employed by Buis and 
Vogt (2J3) will be used here. Equation (136) can be written 
as follows: 



IV < - e *z 2 

dx - e L f 



f W 

a . [ / (a 1 ip 2 + a 2 ^* 2 + a 3 0 2 )dy] • 

[ -W 



W 



1- 



/ [ (l+Ak o )a 4 ^4)-a 4 y0#iJ;-(l-&)Y 0 ^ 2 + (l+Ak o )a 5 e^-a 5 Y0 2 ^]dy 



-W 



w 



-w 



ip 2 + a 2 + a ^ 0 2 ) dy 

Equation (137) can be expressed concisely as 

dV 



(137) 



dx 



£ - R • P • Q 



(138) 



where R = e 2 E^ 2 a is essentially a positive constant. 



r w 

= / ( a ! 'P 2 + a 2 ft, 

J -W 



+ 9 2 ) is greater than zero 



whenever a^ a^ and a^ are greater than zero. It is already 
known that and a^ are positive constants. 



Thus 

zero in order to have P > 0. 



a^ = 1 - (1-3) (1+Ak ) has to be greater than 



63 



From this condition, the first bound for stability is 
obtained 



Ak < 

o 



3 

1-3 



( 139 ) 



in agreement with the linear theory. To prove that Q, 
second expression in brackets in Equation (137), is positive 
definite, the Bun iakows ky- S chwa r z inequality will repeatedly 
be used 

/ fg dy < ( J |f| 2 dy) 1//2 ( f |g| 2 dy) ly/2 (140) 

-a -a -a 

Due to the positive Ak^ insertion, the system’s state 
variables ip , ^ , 9 are positive for any t > 0; therefore, 

the absolute value in the inequality is not necessary. It 
is also known that 

W 



c . 2 ,2 

J Al * 

w 


dy 


< 


Id I 2 


(141) 


. u / A 2'V 

w 


dy 


< 


II 3 I 2 


(142) 


r 2 2 

„/ A 3 0 


dy 


< 


lid | 2 


(143) 



-W' 



State variables in engineering have the properties of func- 
tions in the Hilbert space, and the norm in this space is 
specifically defined as 

a 

( f y 2 dx) 1/2 = ||y|| L (144) 

-a 

Due to the fact that the state variables in the present 
work, ^ and 9, are continuous in space and time domains 

together with their derivatives to arbitrary kth order, it 



64 



is expected that they also have the properties of functions 
in the Banach space. The norm in the Banach space (L^) can 
be defined as: 



( / y m dx) 1/m - 



/ 



= 1 y II, 



for m > 1 



(144a) 



-a 



m 



Then, repeated applications of these inequalities yields 
(1+Ak ) a. 



2wyA 



(t 2 - + 7^)11511- 



A 2 A 3 



1 - 



, CkiB) , _5 

A. 'A 0 A, A,A, . 2 
123 13 



)|! d 



II 3 



min ( 



.2 * . 2’ . 2 
A 1 A 2 A 3 



)||d 



- ,12 



> 0 
(145) 



where the minimum is taken in order to assure the inequality. 
Therefore, for Q > 0, it is necessary that 




> 



(1+Ak ) 
o 



A 

1 




+ 




min 




2WX 



th + d-g) + a s 

A 2 A 3 A 1 A 3 A 2 



) 



(146) 



This is the additional stability condition prescribed by the 
nonlinear theory. The inequality (146) gives the second 
bound for asymptotic stability of the system, which repre- 
sents a spherical surface in the first octant, i.e. ip , 
and 0 >_ 0 . Domains of stability are represented in Figure 
13. The solutions * for *£>(y,T) and 0(y,x) can be obtained 
from Equations (46) and (.57), as follows: 

(y>T) = J* H 1 (t-t') iKyjT 1 ) dx' ( 147 ) 

o 
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w 




1/2 



Figure 13. Schematic of the spherical surface determin- 
ing stability domains for the distributed parameter 
reactor system after a Ak Q insertion. 
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where 



(148) 



A (t-t ' ) 
vZ 



H ^ (x — x * ) = vZ f 3 e 



and 

0(y,x) = j h 2 (t-t ! ) ^ (y , t 1 ) dx 

where _ h (t-t 1 ) 

eZ f CvZ 

H (t-t') = — e 3 

CvE 



(149) 



(150) 



Application of the mean value theorem, in the time domain, 
yields 

[T 

(151) 



i re iqs 

-^(y,T) = ip(y ,t) / H 1 (x-T , )dT' = ifj (y , x) T ± (x) 



where 



BVvE E 

V T > ^ 



, , vE 

(1 - e a ) 



(152) 



and 0 < x < x. 



Also : 



6(y,x) = ijXy.x) / H 2 (x-x') dx'= ijj(y , x) T 2 (x) (153) 



where 



eE 



T,(x) = — r~ (1 - e CvZ a ) 



(154) 



and 0 < x < x . 

Substituting Equations (151) and (153) in Equation (118) 



5II 2 



results in 

~w *vj 

A 1 2 ^ 2 (y,x)dy + / A 2 2 ^2(y, : F)T 1 2 (x)dy + 

y -w -w 

f w 

+ / A 3 2 i|j 2 (y , x ) T 2 2 (x) dy 



(155) 
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Applying the mean value theorem, in the space domain, yields 
|id|| 2 = A 1 2 i(< 2 (y,T) + 2WA 2 2 ^ 2 (y,i) T^x) 

+ 2WA 3 2 !i; 2 (y, : F) T 2 2 (t) (156) 

or 

|i dj| = ^ (y , T ) v / A 1 2 + 2WA 2 2 T 1 2 (t)+2WA 3 2 T 2 2 (t) (157) 

where 0 < T < T and -W < y < W. 

Then the mean value flux can be expressed explicitly as 
f ollows : 

= — ■ ^ -■ (158) 

✓A^ + 2WA 2 2 T 1 2 (t) + 2WA 3 2 T 2 2 (t) 

Typical values for the fast reactor nuclear parameters are 
listed in Table 2, after Solomon and Kastenberg (2j5) . Using 
these parameters, the constants for the system are found to 
be: 



a 2 = 


3.292 


X 


io 13 


2 -2 
cm -sec 




a 3 = 


6.185 


X 


10 27 


-2 0 -l -1 

cm - C -sec 




a 4 = 


9 . 050 


X 


io 9 


-1 

cm-sec 




a 5 = 


1.470- 


X 


io 16 


-2 o-l -1 

cm - C -sec 






2 . 932 


X 


in -26 2 -2 
1U cal -cm 






2 .780 


X 


io' 8 


i 2 - 2 

cal -sec 






1.920 


X 


io 2 


,2 -6 o-l -2 

cal -cm - C -sec 




w = 


3.0 










T 1 (x)= 


70.05 


(l-e“° 


. 0000314x^ 


(159) 


T 2 (t)= 


2.38 : 


K 


io" 12 


-0.0000448X. 
(1-e ) 


(160) 



68 



From Equation (139), the first bound is given by: 



Ak < 1.00331$ 

o 



From Equation (146), the norm, second bound, is obtained 
to be : 



for a Ak^ = 0.156$. The mean value flux from Equation (158) 
is shown in Figure 14. The physical meaning of this mean 



flux provides an order of magnitude information of the in- 
crease in neutron population at the time the surface of 
stability is reached. From Figure 14, it is seen that a 
fast rising flux will quickly reach the surface, whereas a 
more slowly-rising flux will take longer to reach the sta- 
bility domain. 

2 . One-Group Delayed Neutron, Space Dependent 
Step Reactivity Input 

This system is a special case of the previously 
analyzed system. It is desired to find an answer to the 
following question: Given a positive reactivity insertion 

Ak^, is it safer to insert this reactivity uniformly across 
the reactor, Figure 16a, or to insert less reactivity in the 
central region, Figure 16b, provided the total reactivity 
insertion is the same? 

Equation (129) can be rewritten as: 



|| d]| > 2 x 10 5 



value flux i|>(y,T) is depicted in Figure 15. The mean value 




W 




2C £ vE h 
f a 



k 




2C £ vE 
f a 



6^ > dy 



(161) 



vE 
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i|>(y , t) 




Figure 14. Mean value flux vs. dimensionless time plot. 



70 




Figure 15. Mean value flux schematic in the slab 
react or * 
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H I 




Figure 16a 




Figure 16b 



Figure 16 . 



Space-dependent 



Ak 

o 



insertion . 
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Region 



Equation (161) can be rewritten as: 

-2W/3 2W/3 W 

£ -e 2 Z f a [ J* {-} dy + J {-} dy + J* {-} dy (162) 
-W -2W/3 2W/3 



Three steps Ak Q are introduced in the following way: 



k = 1 + Ak 


- y0 


for 


- W < y 


, 2 W 

- 3 


(163) 


k - 1 + Ak 2 


- y0 


for 


2W 

~ - y 


i ¥■ 


(164) 


k = 1 + Ak^ 


- y0 


for 


2W 

— < y 


< w 


(165) 


Substituting Equations 


(163) , 


(164) and 


(165) in 


Equa t ion 



(163) yields, after some algebraic manipulations and intro- 
duction of the constants a,, a„ , a„, a, and a r : 

1 2 3 4 5 



dV 2, 2 

— < - e Z . a 
dx — f 



% 2 + a 2 ^ 2 + a 3 0 2 + (l-3)y9ij) 2 - 



J W/3 
■W , 

- ( 1+Ak j ) ip+ a^yd'&ip- (l+Ak^Ja^S^ + a^y© ^]dy+ 



2W/3 

+ J' [a 3 "ip 2 + a 2 ® )Y 0^ 2 “ ( 1+Ak 2 ) a + 

-2W/3 

+ a^y Q ^ ip - (1+Ak 2 ) a^©^ + a^y© 2 ^] dy + 

+ J [ a 1 1 ^ 2 + a 2 ^ 2 + a 3 0 2 + ( 1-6 ) y Q\p 2 - (1+Ak^ £ ip 

2W/3 . 

4 a^yG^ijj -(1+Ak^)a^ 0 4* a^ y0 2 ^]dy 



4 

(166) 



where 



a 1 1 = 1 - (1-6)(1+Ak 1 ) (167) 

a 1 " = 1 - (1-6)(1+Ak 2 ) (168) 
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Following the same steps as in the previous analysis, the 
following conditions for asymptotic stability are obtained: 



Ak 



1 



< 



3 

1-6 



Ak 2 



< 



6 

1-3 



conditions which agree with the linear theory, and 




(1+Ak ) a, 

— (— + 

A 1 A 2 



5 

-7—) - min ( 

A 3 






wl ( !l. + 

3 A 2 A 3 



(1-3) 

A 1 A 3 




(169a) 

(169b) 



(170) 



3 2 a 3 



d||. 



(1+Ak») a, a_ a 

A ^1T + A^ } " m±n (_J T’ 2* 2 

A 1 A 2 A 3 A, A„ A, 



) 



4W 



a _4 ( 1-3) a 5 . 

3 A, S A„A_ A, A, .2' 

1 2 3 13 A ^ 



(171) 
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These are the additional condition of stability required by 
the nonlinear theory. In order for the total reactivity 
insertion to remain the same for both strategies, one must 



have 



-W 



-2W/3 



2W/3 



2 « J 



Ak dy = 
o J W 



W 



/ Ak i dy + h / 

-W Ot.T / 



Ak 2 dy 



(172) 



-2W/3 



From Figure 17, it is seen that in Region 2 in which Ak 2 < 
Ak Q , ||d|| 2 > || d || , so that the region of stability is attained 
sooner than in the case of uniform Ak Q insertion; but ||d||^> 
||d||, so in Region 1 of the reactor, where Ak^ > Ak^, the 
region of stability is attained later than for the case of 



uniform Ak insertion. However, the flux in the central 
o 
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Figure 17. Spherical surfaces for the case of space- 
dependent Ak^ insertion. 
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region of the reactor plays a more important role in the 

transient behavior of the reactor; consequently, the space 

dependent reactivity insertion, with Ak^ < Ak^, is safer 

than the uniform Ak insertion. 

o 

3 . One-Group Delayed Neutron, General Reactivity Input 
For this case, Equations (37), (46), (57) and (66) 

are the governing equations. The feedback coefficients, y, 
and y , are assumed to be negatives; and y^ are the 

Doppler and Expansion coefficients of reactivity, respective- 
ly. The same Liapunov function, Equation (122), is used 
and a similar development previously discussed leads to the 
following equation, in analogy to Equation (136): 



dV 

dx £ -e 2 E f 2 a 









[ a ^ 2 + + a 3 e 2 + (1-3H 2 F(0) + 

+ a 4 l|)F(0) a 5 0 iJjF (0 ) - ( l+Ak Q ) a ^ \p- ( l+Ak Q ) a 5 0tf> ] dyj ( 173 ) 
where the feedback term is 

f(9) ■= It 1 0 + -jrj— It 1 " 11 



- t 1_n ]+ 

o 



IyJp 



f _ [T l +m _ T0 l +m](174) 



here R and P are constants equal to 300°K. F(0) is positive 

for all values of n > 0 and m > 0. It is noted that the 
Doppler feedback is the primary feedback mechanism in large, 
ceramic fueled fast reactors, and expansion is usually the 
dominant feedback in small metal fueled fast reactors. For 
n=l, Equation (66a) is used. A different approach, more con- 
servative, is used to analyze Equation (173). The expression 
inside the curly brackets in Equation (173) has to be posi- 
tive definite in order to ensure asymptotic stability. With 
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this in mind, a comparison of the orders of magnitude among 
the terms is carried out. 



d V 
dT 



< 



e 2 Z 



a [a^ 2 +^{a 2 ^ 



( 1+Ak^ ) a } + 



+ 0 {a 3 0 - (l+Ak o )a 5 i^) + iJjF (6 ) { (1-3) ip+a^f+ a 5 0}]dyj- (175) 

In order to get a positive definite function inside the 
curly brackets, the following conditions have to be satis- 
fied: 

o 

(a) a- > 0 which yields Ak < ■= — 5 - 

1 ' o 1-p 

(b) a 2 'fo “ (l+Ak Q ) a^ijj > 0 which yields 

4 < 

f> ^ 1+Ak Q ) a 4 

and 



(c) a^0 - ( 1 +Ak^) a<- ip > 0 which yields 

± < Ji 

e (l+Ak )a, 

O 3 



Typical values for the feedback coefficients are: 
y= 10 _ 4 /°K, y= 7xlO - 6 /°K, y= 9.5 x 10 _ 6 /°K. Evaluating 

i) £j 

these conditions for a typical fast reactor, with a Ak^ = 



156$, 


it 


is 


obtained that: 


(a) 


Ak 


o 


< 1.00331$ 


(b) 


t 


< 


3634.0 




c 






(c) 


i 

e 


< 


4.203 x 10 11 



These three conditions have to be satisfied simultaneously. 
The domain of stability is shown in Figures 18a, 18b and 

18c. It can be seen that the domain of stability is very 
large when the coefficients of reactivity are negative. 
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Figure 18a. Stability domains for the distributed 
parameter reactor system after a general reactivity 
insertion. 
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Figure 18b, Stability domains for the distributed 
parameter reactor system after a general reactivity 
insertion . 
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Figure 18c, Stability domains for the distributed 
parameter reactor system after a general reactivity 
ins er t ion • 



80 



IV, CONCLUSIONS 



Liapunov's Second Method has proven to be a useful 
method to find stability domains in Nuclear Reactor Control, 
The fact that a Liapunov Function is not unique makes the 
choice of such a function to depend upon the experience and 
ingenuity of the analyst. Once an appropriate choice of the 
Liapunov Function has been made practical results can be 
found . 

For the Point-Kinetics Model, a Liapunov Function was 
obtained without imposing an a priori bound, as most of 
the conventional methods do, Domains of stability are pre- 
sented for different cases using typical data from thermal 
reactors . 

The Distributed-Parameter Reactor System was studied 
using the concept of a norm, an a priori bound, and the 
variation with time of this norm, defined in a Hilbert 
space, was analyzed to obtain a domain of stability. Data 
from a typical fast reactor was used to define the spherical 
surface of stability and to estimate a mean value flux, in 
space and time, when this surface is reached. The case of 
space dependent Ak^ insertion, with less reactivity in the 
central region, was found to be safer than the case of 
uniform Ak^ input. 

The linear theory stability condition appeared through- 
out the analyses giving confidence to the obtained results, 
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in addition to that condition and additional condition(s) 
appeared in each case, required by the nonlinear theory. 

The inclusion of the multigroup delayed neutrons is 
recommended for further improvement of this analysis. 
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APPENDIX A 



THE EIGENVALUE INEQUALITY 



Let X > 0 be the least eigenvalue of the system 



ill 

dy 2 



+ Xijj = 0 



(Al) 



with boundary conditions: (a) ^(+W) = 0 



the following inequality can be stated: 

w / W 

\ / i^dy < / (f^) 2 dy ( A2 ) 

Proof of this inequality can be done using the calculus 
of variations . 

Let 



W 

J XiJj 2 )dy 0 






Applying Euler’s equation: F^,= 



(A3) 

(A4) 



Equation (Al) is obtained. 

Solving the eigenvalue problem X is obtained to be: 

(A5) 



X = — = b 2 l 2 

4W 



Then 



-W 



(^) 2 dy > B 2 L 2 






^ 2 dy 



(A6) 



-tv 
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